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Abstract 


In order to find the attitude of a spacecraft with respect to a reference coordinate system, vector 
measurements are taken. The vectors are pairs of measurements of the same generalized vector, taken in 
the spacecraft body coordinates, as well as in the reference coordinate system. We are interested in 
finding the best estimate of the transformation between these coordinate systems. The algorithm called 
QUEST yields that estimate where attitude is expressed by a quaternion. QUEST is an efficient algorithm 
which provides a least squares fit of the quaternion of rotation to the vector measurements. QUEST, 
however, is a single time point (single frame) batch algorithm, thus measurements that were taken at 
previous time points are discarded. 

The algorithm presented in this work provides a recursive routine which considers all past 
measurements. The algorithm is based on the fact that the, so called, K matrix, one of whose 
eigenvectors is the sought quaternion, is linearly related to the measured pairs, and on the ability to 
propagate K. The extraction of the appropriate eigenvector is done according to the classical QUEST 
algorithm. This stage, however, can be eliminated, and the computation simplified, if a standard 
eigenvalue-eigenvector solver algorithm is used. The development of the recursive algorithm is 
presented and illustrated via a numerical example. 

L INTRODUCTION 


The problem of finding attitude from vector observations is stated as follows. A sequence, b., i=l,2, 

... , k of unit vectors is given. These unit vectors are the result of measurements performed in 
vehicle cartesian coordinates, of the directions to known objects. The sequence, r., i=l,2, ... , k of 

unit vectors, is the sequence of the corresponding unit vectors, resolved in a reference cartesian 
coordinate system. We wish to find the attitude matrix, A, which transforms vectors from ^he reference 
to the body coordinates. Obviously, A has to be an orthogonal matrix. In 1965, Wahba 1 posed this 
problem as a least squares problem as follows. Let 


WA> *ii5jlV ^il 


2 


( 1 ) 


find that orthogonal 3x3 matrix, A, that minimizes L. We can weigh each measurement separately 
according to the accuracy of the particular vector measurement. In addition, we may want to find the 
quaternion, rather than the matrix, representation of attitude. In such case (1) is replaced by 
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( 2 ) 


k 

J(q) = \ ^Jb.- A(q)r.| 
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where a., i=l,2, ... , k are the positive weights assigned to each measurement. In (2) we are looking 

for that quaternion, q, which minimizes J. Note that instead of minimizing J, we can maximize g defined 
as 

g(q) = 1 - J(q) (3) 

It can be shown 2,3 that g(q) can be written as 

g(q) = q T K q (4) 


where K is constructed as follows. Define 


m. 


= I a 
i =1 i 


(5.a) 


Then 


I t 

a = — I ab r 
m i = 1 i i i 

k 

1 k T 
B = — X a.b.r 
m 1=1 i i i 
k 


S = B + B 
1 k 

z = -J- ,X a.(b. x r.) 
m i=l i i i 
k 


(5.b) 


(5.c) 

(5.d) 

(5.e) 


' S-oI 

z 

T 


Z 

a 


( 6 ) 


2 3 . ... 

where I is the third order identity matrix. It was shown ’ that q*, of unity length, which maximizes 
g(q) in (4), satisfies the equation 


K q* = X q* (7) 

where X is a, yet undetermined, Lagrange multiplier. We realize that X is an eigenvalue of K and q* is 
the eigenvector which corresponds to X. Substitution of this solution into (4) yields 

g(q*) = X (8) 

and since we wish to maximize g, we choose X , the largest eigenvalue of K, as the desired 

max _ 

2 

eigenvalue, and then, q* is the eigenvector which corresponds to this Davenport showed that once 

X is found, there is no need to solve for the eigenvector of K, since y\ the optimal vector of 
4 5 

Rodrigues parameters , (also known as Gibbs vector ) can be computed as follows 
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( 9 ) 


r = [(X + o)i -s^z 

max 

and the optimal quaternion can be found using the known relation 


M 



Shuster ’ showed how to, easily, compute X to arbitrary accuracy, and how to deal with a singular 

matrix in (9). It was also shown there that X is close to 1 and is exactly 1 when the measurements 

max 

are error-free. (This property is due to the fact that all a.’s in (2) add up to 1, or, equivalently, 

the introduction of the normalizing factor m in to (5.b, c, and e)). The algorithm for obtaining X 

* max 

and q* from vector observations discussed above is known as the QUEST algorithm. 

QUEST is a single point attitude determination algorithm; that is, it utilizes the vector measurements 

obtained at a single time point and uses them, and them only, to determine the attitude at that time 

point. This way, the information contained in past measurements, is lost. This fact has been recognized 

and in 1989, Shuster presented an algorithm which he named Filter QUEST that processes vector 
measurements recursively. The Attitude Profile Matrix, B, defined in (5.c), which plays a central role 
in the algorithm, is updated recursively for use in the QUEST algorithm. Much attention is given, in 
that paper, to covariance calculations. 

In the present work, the matrix, which is updated recursively, is the K matrix defined in (5) and (6). 
Indeed, as can be seen in the algorithm described above, K is the most important element in QUEST. In 
the following section, we start our presentation of REQUEST by considering, first, the recursive 
time-invariant algorithm. Then, in Section III, we develop the recursive algorithm for the time-varying 
case and present an example. In Section IV we list the algorithm in a unified form. Finally, in Section 
V, we present our conclusions and recommendation for further work. 

IL THE RECURSIVE TIME-INVARIANT ALGORITHM 

Assume that the body axes are non-rotating with respect to the reference axes. Also assume that k 
vectors have been processed using the QUEST algorithm. 

Let 


k 

m = la 
k i = l i 


where m^ is not necessarily equal to 1. Also define 


(11 -a) 


1 T 

o = — L ab r 
k m, i=l i i i 
k 


1 T 

B, = — I abr 
k m, i = l i i i 
k 




(ll.c) 

(II.d) 
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(He) 


1 

z = — £ a.(b.x r) 

k m k i = 1 1 1 J 


The parameters o k> S k and z^, are then used to compute K, which for the case of k measurements, is 

denoted by K . The latter is used to find the optimal quaternion, based on k pairs of measured vectors. 
* k 

(Note that QUEST doesn’t require the computation of itself). As mentioned earlier, the coefficient 

m is used in (ll.b, c, and e) to normalize the weights, a., such that X^ is closed to 1. (See Ref. 5 

for the solution of X ). Now suppose that an additional measurement has been acquired; that is, the 
max 

k+lsr pair has to be processed. The question is then, do we have to re-compute the K k+J matrix anew or 
can we, perhaps, update K k to included the added pair. As will be shown next, the latter is possible. 

In fact, it forms the basis for the REQUEST algorithm. We formulate this quality of K in the following 
proposition. 


Proposition 1; Let 


^*k+l a k+l^k+i r k+l 


5B = a b r 
k+l k+1 k+1 k+1 


5S = 8B + 8B* f 
k+1 k+1 k+1 


(12.i) 

(12.b) 

(12x) 


5z k + l = “k+A+i* r k+i> 


(12.d) 


5K = 
k+1 


8s k .r 'Vi 1 

8z k+i ' 

8z k 

L k + 1 

8o k+i - 


(12.e) 


that 


m k 1 

K * — K + — — 6K 

k+1 m km. k+1 


k+l 


k+1 


(13) 


Proof: By definition 


Now, it can be easily verifies that 


k+l 

m = £ a. 
k+l i=i i 


(14.a) 


a i k+l t 1 k T 

c = — — I a.b T r. - — — X a.b'r. + 
k+i m i = i iii m i«i i i i 

k+l k+l 

. , k +l _ , k 

B = E.a.b.r = I.a.b.r! + 


1 .T 

a b r 

m k+l k+l k+1 
k+l 

1 


. a b r 

k+l m i = l“i”i i m i = l i i 1 m . k+1 k+1 k+1 


k+1 


k+1 


k+1 


(14.b) 


(14.c) 
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(14.d) 


s u , = B u , + B I , = K + B l + SB + 8B T 
k+1 k+1 k+1 k k k+1 k+1 


k+1 


k+1 m i = l i' i i' m 

k+1 k+1 


A 1 1 k 1 

z - - ~ .Z,a.(b.x r.) = — .£a.(b.x r.) + a, (b x r ) 

j = i i i r m. . k+i k+i k+i 


(14.e) 


k+1 


Using the definitions in (11) and (12), (14.b) to (14.e) can be written as follows 


t - m fc „ 1 - 

k+i m. . k + m. . k+i 


k+1 

m. 


k+1 

1 


8S. 


k+l m, , k m. k+i 
k+1 k+1 


(15.a) 


(15.b) 


m 

k 1 - 

z = z + 8 z 

k+l m, k m, k+l 
k+1 k+l 


(15.c) 


Whe “ a k+l’ S k+l’ and z k+i def " med “ ( 15 ) m used to form K k+j , using the format of (6), (13) results. 
This ends the proof. 


We have assumed here that we add only one new measurement to the k measurements that were already 
processed. This can be extended to the case where two or more measurements are added as a group of 
measurements. Suppose that K was computed n times where at each time, one or more measurements were use 
to compute (initially) or to update K. Let this K be denoted by n=l,2, ... , where is computed 

using (5) and (6), and where the index k is the number of measurements used to compute K . Suppose that 

j new pairs of vectors are measured and we want to use them in the updating of K. We can, of course, 
update K n j times, using the algorithm presented in Proposition I, and obtain the updated K, or we can 

lump the new j measurements together, and update K only once. The latter is performed according to the 
algorithm listed in the following proposition. 


Proposition 2: Let 


8m 

n+1 


k*j 

i-f.l *1 


(16.a) 


where k is the number of, already processed, pairs of vector measurements. 


k+j T 


8a = . E a.b r 

n +1 i=k+l i i i 

(16.b) 

k + j T 


8B = E abr 

n +1 i=k+l i i i 

(16.c) 

8S = 8B + 8B T 

n +1 n +1 n +1 

(16.d) 

k+j 


8z = . E a.(b x r ) 

n +1 i=k+l l i i 

(16.e) 
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(16.0 


8K = 

n+l 


5S -5a I 

n+l n+l 

8z 

n+l 

5z T 

L n + l 

5a , 

n+l J 


then 

and 


m 


n+l 

m 


K 


n+l 


m 


n +1 


= m + 5m , 

n n+l 

(17.a) 

K + 1 8K 

n m n+l 

n + l 

(17 .b) 


This proposition can be easily proven along the lines of the proof of the first proposition. The case 
described in the first proposition is a special case of the latter, we chose to split the introduction 
of the updating of K into two cases merely for methodical reasons. 


HI. THE RECURSIVE TIME-VARYING ALGORITHM 

The updating algorithm of the static case can now be extended to the case where the body rotates 
between measurements. In the ensuing development we distinguish between two cases; namely, the 
error-free propagation case, and the propagation which is based on angular rate measurement, and as 
such, is contaminated by rate-measurement errors. 


ni l Error-Free Propagation 


Assume that at time t , k pairs were processed, then the body rotated to a new orientation and there, 

at time t , j new vector measurements were performed. We wish to find the least squares fit of the 

quaternion to the first k measurements, at this new time point, and then do the same when the new j 
measurements are considered too. So first we are interested in finding q* +1M which is the quaternion 

that expresses best the attitude at time t^, based on the first k measurements which were performed 

previously, at time t . Let us re-write the cost function of (4) for q at time t^ based on the first k 

measurements which, as mentioned, were performed at time t^ 

g(q ) = q T , K , q „ ( lg ) 

K '' M n/n M n/n n/iiTi/n 

It is well known 9 that during the rotation, q changes according to the differential equation 

q = i0q <•«> 

where Q is a 4x4 skew symmetric matrix whose elements are the body components of the vector of the 
angular velocity of the body with respect to the reference frame. The solution of (19) yields 


q(t .) = ®(t ,t )q(t ) 

n+l n+l n n 


( 20 ) 


Ideally, when £1 is known perfectly, the matrix ®( t n+1 * t n ). known as Transition Matrix, transforms the 
quaternion which represent the attitude at time t , to that which represents attitude at time t^. For 

simplicity of notations, we denote it, simply, by O. The quaternion which we wish to transform from 

time t to time t is q , thus we set q(t ) = q . . Finally, we denote the quaternion, to which 
n n+l nM n n/n 

q is transformed, by q , , , thus we set q(t ) = q . Consequently (18) becomes 
H n /n n+l/n n+l n+l/n 


q = q , 
M n+l/n nto 


( 21 ) 
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( 22 ) 


Since £1 is skew symmetric, O is orthogonal, thus we can write 




1 _ T 

q = O q 
n+l/n 11 +IA 1 


Substitution of q of (22) into (18) yields 

n/n 


“w = - <„/» ® K.,y v,„ < 23 > 

We realize that the problem of finding q that maximizes g has been transformed into the problem of 

n/n r 

finding q which maximizes g’. Let 
n+lAi 


then (23) becomes 


K 


n+1 /n 


<E> K 0 T 

n/n 


(24) 


^ ^n+l/n^ ^n+l/n ^n+l/n ^n+l^i 


(25) 


One may ask oneself whether the problem of finding q , which maximizes g\ is still related to 

n+l /n 

Wahba’s problem; that is, will the maximization of g’ yield a quaternion which is a least squares fit 
to the k vector measurements. The answer is, of course, positive, since the maximizing q is 

n+l /n 

directly related, through (22), to q which maximizes (18), and the latter is the solution of Wahba’s 

n/n 

problem, given the k measurements. It can be shown (see the Appendix) that, like before, q* , which 

n+l/n 

is the q n+j ^ that maximizes g’, given in (25), satisfies the equation 


K q* = X q* C261 

n+l/n M n+1/1 n+l/n M n+l/n ^ 

and that q* is the eigenvalue of K which corresponds to the largest eigenvalue of K .It 
n+i/n n+l/n ** n+l/n 

is interesting to note that this solution to the constrained optimization problem is not specific to 

attitude determination. It stems from the fact that the cost function is defined as a quadratic form of 

a square matrix and that q is required to be of unity length. (See the Appendix). Also note that 

although we assume error-free propagation, the measured vectors contain measurement errors. Finally, 

note that K n+1/n » being a result of a similarity transformation on K , has the eigenvalues of the 

latter even though its eigenvectors are different. 


Now that we have established the fact that K n+j/n is the proper K matrix for finding the least squares 
fit of the quaternion at time t based on all past k measurements, we can include j more measurements 
performed at t^. For this we use (17.b) of Proposition II. Consequently from (24) and (17.b) we 
obtain 


K " =OK <p 

n+l/n n/n 


(27.a) 


m . 

K = — K + — — SK 

n+l /n+l m . n+l/n m n+l 

n+l n+l 

We demonstrate the algorithm by way of the following example. 


(27.b) 
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Example: 


Data base : 

Given ar$ 4 error free vectors in the reference coordinate frame: 

0.267 "I [-0.667 1 [ 0.267 1 -0.447 

0J35 r2 = -0.667 r3 = -0.802 r4 = 0.894 

0.802 -0.333 L 0535 J L 0 000 

and a rotation from the reference to body axes described by the following Euler vector: 



((> T = [0.9, 0.2, 0.8] 

The corresponding quaternion is: 

q(l) T = [0.423, 0.094, 0.376, 0.819] 

The four r vectors are transformed to the body frame and noise is added to the transformed vectors. The 

noise elements added to each component of the transformed vectors is drawn from a random number 

generator. The standard deviation of the noises are: 

al = 0.01 a2 = 0.05 o3 = 0.03 a4 = 0.02 

The noise element added to each component of r. is drawn from the random number generator whose 

standard deviation is a., i = 1, 2, 3, 4. The vectors are then normalized. The resulting simulated 

measured vectors in body frame are then: 

-0.280 ' 

b3 = -0.030 b4 
0.959 


'-0.985 
b2 = -0.120 
-0.123 



and the weights are chosen to be: 

a = a.' 2 i = 1, 2, 3, 4 

i i 

Application of QUEST to the first two pairs : 

Using, initially, at time t , the first two pairs of vectors, rl and bl, and r2 and b2, we obtain K^. 
Its largest eigenvalue and the corresponding eigenvector, which, according to our notations, is q^, 
are: 

X = 1.0003551 = [0.427, 0.105, 0.383. 0.813] 

The corresponding transformation matrix, A , the correct matrix, A(l), and the difference (error) 
matrix, are: 
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'0.685 0.712 0.156 ' 


"0.700 0.695 0.164' 

A ,,= 

-0.532 0.343 0.774 

A(l) = 

-0.536 0.361 0.763 


0.497 -0.613 0.614 


0.471 -0.622 0.625 


A - A(l) = 

l/l v ’ 


-0.015 0.017 -0.007 
0.004 -0.018 0.011 
0.026 0.009 -0.012 


The Euclidean norm of the error matrix is: 


|A W - A(l) | = 0.044 

This error steins, of course, from the measurement error contained in the b vectors. 
Rotation of the body coordinate system: 


(28) 


We assume that after processing the first two pairs, which yielded A^, the body rotates for 1 sec at 
the following angular rate: 

co T = [0.1, 0.2, -0.3] rad/sec 

The matrix <l> which propagates the quaternion of this rotation (see (21)), and, A A, the attitude matrix 
which expresses the change in the body coordinates are: 


<D = 


0.983 -0.149 -0.099 0.050 
0.149 0.983 0.050 0.099 
0.099 -0.050 0.983 -0.149 
-0.050 -0.099 0.149 0.983 


AA = 


0.936 -0.283 - 0.210 
0.303 0.951 0.068 
0.181 -0.127 0.975 


Measurement update of K: 


We use AA to transform and to the new time point, t^. Using these b*s as the simulated 
measurements at t 2> we compute according to (16), and update K, using (27), as follows: 


K = 3> K O 1 
2/1 1/1 

m 

K 2/ 2 = k 2 /i + ^ 5k 2 


The largest eigenvalue and the corresponding eigenvector (which is of are: 


X V2 = L0001957 \ a = [°- 4 °2. 0.253, 0.282, 0.834] 
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The corresponding attitude matrix is: 


f 0.713 0.673 -0.195 



-0.267 0.518 0.813 
0.648 -0.528 0.549 


Check : 

We wish, now, to check this result. This is done as follows. We use AA to transform and to the 

new time point, (Recall that b 3 and b 4 were already transformed in order to compute 8K 2 ). Now we 

apply the QUEST algorithm to all four pairs of r and b. The resulting quaternion should be equ^l to the 
quaternion updated by the REQUEST algorithm. Indeed the two quaternions agree to at least 10 . 

Remark: 


When we compare, A(2) = AA-A(l), which is the correct matrix which transforms from the reference to 

body axes at t , to the attitude matrix A , obtained by REQUEST (and, as we just checked, by QUEST as 
2 */ ^ 


well), we obtain: 


0.005 - 0.006 0.000 



A(2) = 


-0.001 0.007 -0.005 


-0.006 -0.001 0.007 


The Euclidean norm of this error matrix is: 

|A W - A(2)1 = 0.015 

This error stems from the measurement noise in the b vectors and not from the algorithm. We note that 
the latter error is smaller than | A - A(l) | shown in (28). This is expected, since A ^ is computed 

using four pairs of vectors whereas A is computed using only two. 


E0L2 Noisy Propagation 


In the preceding developments we considered the presence of noise only in the measurements and assumed 
that the angular rate vector, to, was known to us perfectly. We wish to consider now errors also in our 
knowledge of to. Let us denote the measured, or computed, to by co^. We also assume that the error is 


additive, thus we can write 


at = co + e 

m 


(29) 


where e is the error component in the measured angular rate vector. We distinguish between two cases; 
namely, short time application, and long time application of REQUEST. The two are treated next. 


TTT.2. 1 Short Mission Duration 

Since a typical update rate is once per second, typical gyro noise does not cause a considerable 
attitude error during such a short period. In fact, even with an update rate of once per 10 seconds, 
the attitude error amounts to a very small attitude error. To illustrate this point, we turn to Jhe 
example. Suppose that we use a triad of single axis gyros, each having a constant drift rate of 1 /h, 
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which is about 100 times larger than that of inertial grade gyros. We use the first three measurements 
to compute and from it, A^. We then propagate using the gyro-error ridden transition 

matrix, and obtain K and then compute the corresponding attitude matrix, A . In parallel we do 

z/1>m 2/1, m 

the same using ®, the correct transition matrix, and obtain A , the corresponding attitude matrix. 

Doing so, we discover that the largest difference between the magnitude of the elements of the two 
attitude matrices, m and A^, is less than 5.23-1C . Next, following the REQUEST algorithm, we 

use the fourth measurement at time t to compute 8K , update both, K and K , and compute the 

corresponding attitude matrix for the correct and erroneous propagations. Hie largest error between the 
elements of the, updated, two attitude matrices is less than 2.55- 10 . We see two interesting facts. 
First, indeed, the gyro error has little effect on the propagated K and, consequently, on the, 
propagated and the updated, attitude matrices. Second, the incorporation of a new measurement reduces 
the little error, caused by gyro drift, even further. As a consequence of this discussion, we conclude 
that for a short mission duration the build up of attitude errors as a result of gyro drift is 
negligible and the algorithm given in (27) is adequate. 


III.2.2 Long Mission Duration 

Space missions where QUEST is traditionally being used, are of long duration, therefore the initial 
measurements are propagated through the repeated use of (27.a) to the current time. This in turn 
reduces the accuracy of those measurements, and as time goes by they may corrupt the attitude rather 
than improve it. Consequently, we wish to gradually reduce the influence of old measurements, and 
eventually eliminate them altogether. This is usually done using the Fading Memory 1 " concept. 
Accordingly, instead of using (27 .b) for updating K, we may want to use the following algorithm 


m 


K = p — 

n+l/n+1 r n m 


K 


n+l/n m 
n+1 n+1 


5K 


n+1 


(30.a) 


where 0 < p n < 1. Note that p n has to be larger than 0 for (28.a) to yield a meaningful K when only 
one measurement is performed at t . Also note that when no process noise is present, we set p = 1 

which keeps the same relative weighting of past and present measurements as in (27.b). The value of p 
can be determined experimentally where a larger propagation noise is compensated by a smaller p value. 
Note that can vary from step to step allowing the consideration of changing gyro noise. It should be 
noted that the introduction of m in the REQUEST algorithm stems from our wish to maintain X si. 

n max 

This is important if we use the classical QUEST method for solving for X [see ref. 7]. (If we use a 

max 

given eigenvalue-eigenvector solver routine, this is irrelevant). When (28.a) is used and we are still 
interested in having X close to 1, we have to replace (28. a) by 


p m 


_ £ ** K -f ^ 

n+1 /n+1 p m + 5m n+l/n p m + 8m 

n n n+1 r n n n+1 


8K 


n+1 


(30.b) 


Note that, as before, this K update algorithm still assures proper weighing of the measurements; that 
is, the measurement noise is properly considered. 


IV. ALGORITHM SUMMARY 


The REQUEST algorithm is summarized as follows. 

1. Use the k measurements performed at the starting point, t jf tocompute K . First compute: 
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(31.a) 


m k = i=i a i 

a = — I a.b T r. 
1 = 1 i 11 

1 k T 

B = — E a.b.r. 

m k 1=1 i i i 


Then compute: 


2. Form the angular rate matrix: 


S = B + B 
1 k 

z = — .E a.(b.x r.) 
m 1=1 i i i 
k 


S-oI 




'1/1 

L ‘ T l 

1 a J 


0 

CO 

-CO 

CO 


z 

y 

X 

-CO 

0 

CO 

( 0 

z 


X 

y 

CO 

-CO 

0 

(0 

y 

X 


z 

-CO 

-CO 

-© 

0 


(31.b) 


(31.c) 

(31.d) 

(31-e) 


(31.f) 


(32) 


where, o), i= 1,2,3 are the components of the body axes, angular rate vector, 
i 

3. Compute 4*, the transition matrix from time t^ to time t^, corresponding to this, generally 

time-varying, angular rate matrix. (Algorithms^ for computing O can be found in standard Control 
Theory or State Estimation texts. See e.g. Gelb .) 


4. Propagate according to: 

K 2 „ - <■> K 1( / 

5. Compute 5K as follows: 

2 k+j 

8m = E a. (34.a) 

2 i=k+l i 

(where k is the number of, already processed, pairs of vector measurements, and j is the number of 
new measurement pairs performed at time t ). 


J? 

II 

II K* 

kM + 
+ «-• 

(34.b) 

k+j T 

8B = E a.b.r 

2 i=k+l i l l 

(34.c) 

5S„ = 8B„ + 8bI 
2 2 2 

(34.d) 

k+j 

8z = X a.(b x r.) 

2 i=k+l l i i 

(34.e) 
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8K 2 = 


5S„ - 6ct I 

Sz 

2 2 

2 

K 

5oJ 


then set in the range 
and compute 

K 

2/2 


0 < p ^ 1 
n 


p m 
K i l 

P i m i + 5m 2 


K + — 5 — 5K 

2/1 p m + om A 2 
11 2 


In preparation for the next time update, compute 


m = m + 8m 
2 1 2 


(34 .f) 


(34.g) 


(34 i) 


(34 .i) 


6. Only if there is an interest in extracting the attitude from K^, compute the attitude at this time 


point (t ), otherwise go to step 7. 


The extraction of attitude from K can be done using the 

2/2 


algorithm given in QUEST, or any standard software package that can compute eigenvalues and 

eigenvectors of a symmetric matrix (e.g. Matlab™ or Mathcad™). If the latter approach is chosen, 
then, first, select the largest eigenvalue of K^, and, then, compute the corresponding 

eigenvector. 


7. Go to step 2 and increase the appropriate indices by 1, or stop if so desired. 


V. CONCLUSIONS AND RECOMMENDATIONS 

In this work we presented a recursive algorithm for attitude determination, from vector observations, 
that was derived from QUEST. The new recursive algorithm, which we call REQUEST, is based on the 
propagation and update of the K matrix, one of whose eigenvectors is the sought a t tit ude quaternion. 
Using REQUEST, we do not lose information gathered by measurements performed at previous time points, 
and since we use prior information, even one measurement at a particular time point, to which K is 
propagated, is sufficient for updating the attitude. We showed how to apply the algori thm to cases 
where more than one measurement is taken at the new time point. We demonstrated that under normal 
conditions, and for short mission durations, there is no need to treat propagation noise (also known as 
process noise). For long mission durations we do have to consider the process noise. This is done using 
the Fading-Memory notion whereby the weight of the contribution of old measurements to K is reduced 
with time. We presented an example to illustrate the algorithm. 

As mentioned, the new algorithm shows how the propagate and update K, but once K is computed, its 
largest eigenvalue and the corresponding eigenvector, which is the sought quaternion, are found using 
the method of QUEST. If, however, a standard eigenvalue-eigenvector solver algorithm, is used, then the 
eigenvalue and eigenvector can be found directly without solving for Rodrigues parameters, and without 
the need to be concerned about matrix singularity problems (see (9)). 

As a follow up to this work, it is recommended that REQUEST be tested using real spacecraft data, and 
be tested against other recursive algorithms, such as the extended Kalman filter. 
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Appendix 


In this appendix we prove that a cost function formulated as a quadratic form of a real symmetric 
matrix, with a unity constraint on the vector part of this form, has the following two qualities: 

I. Its maximum is equal to the value of the matrix largest eigenvalue. 

II. The vector which maximizes the cost function is the matrix eigenvector which corresponds to this 
eigenvalue 

We present the proof in a form of a question and an answer as follows. 


Problem: Given 

p = x T M x (A.1) where |x| = 1 (A.2) 

and M is an nxn symmetric matrix, find x which maximizes (J, 

Solution: We use the method of Lagrange multipliers to incorporate the constraint of (A.2) in the cost 
function expressed in (A.1). Accordingly, we wish to maximize <p(x) given by 

cp(x) = X T M x + A(1 - x T x) (A.3) 

We denote the maximizing x by x*, then we can express x as follows 

X as x* + eh (A,4) 

where € is a scalar. Substitution of the latter into (A.3) yields 

<|>(e) = (x* + eh) T M (x* + eh) + Ml - (x* + eh) T (x* + eh)] (A.5) 

An extremal point of 4>(e) satisfies the following 

d< ^| =0 for all h (A.6) 

de 1 e=0 

Now it can be easily verified that since M is symmetric, 

^1 _ = 2h T (Mx* - Ax*) (A.7) 

de 1 e=0 

Application of the condition for a stationary point of (A.6) to (A.7), yields 

h T (Mx* - Ax*) = 0 for all h (A.8) 

The latter condition can be met if and only if 

Mx* = Ax* (A.9) 


Substitution of Mx* given by the last equation into (A.1) yields 

ll = Ax* T x* 

max 

and since x* is of unit length, x* T x* = 1, therefore (A.10) becomes 

Li = A, 

max 

and u takes its maximal value when A, is A , which is the largest eigenvalue of 
“max max 

M. (Note that since M is symmetric, its eigenvalues are always real). Then 

H = A, 

max max 


(A.10) 

(A.11) 


(A.12) 


and x* is the eigenvector of M which corresponds to A 

® max 
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